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ABSTRACT 

A probabilistic evaluation of an eight-ply graph i te/epoxy quasi-isotropic 
laminate was completed using the Integrated Composite Analyzer (ICAN) in con- 
junction with Monte Carlo simulation and Fast Probability Integration (FPI) 
techniques. Probabilistic input included fiber and matrix properties, fiber 
misalignment, fiber volume ratio, void volume ratio, ply thickness and ply 
layup angle. Cumulative distribution functions (CDF's) for select laminate 
properties are given. To reduce number of simulations, a Fast Probability 
Integration (FPI) technique was used to generate CDF's for the select proper- 
ties in the absence of fiber misalignment. These CDF's were compared to a 
second Monte Carlo simulation done without fiber misalignment effects. It is 
found that FPI requires a substantially less number of simulations to obtain 
the cumulative distribution functions as opposed to Monte Carlo Simulation 
techniques. Furthermore, FPI provides valuable information regarding the sen- 
sitivities of composite properties to the constituent properties, fiber volume 
ratio and void volume ratio. 


SYMBOLS 


E cxx 

E fll> E f22* G fl2 

Em^m 

E (x),Pf 

f(x) 


fx(X) 

G cxy 

g 


composite elastic modulus (Mpsi) about structural axes 

fiber elastic moduli about material axes 

matrix elastic moduli 

cumulative distribution function 

probability density function 

joint probability density function 

composite shear modulus (Mpsi) about structural axes 

FPI limit state function 
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kf ,k m ,k v 


u 

x 

y 

z 

a,0 

a cxx 

r(x) 

X,ic 

lj,a 

v f 12 » v m 


linear approximation of g 

incomplete quadratic approximation of g 

fiber volume fraction, matrix volume fraction, void volume 
fraction 

standardized normal deviate 
uniform deviate 

normal, Weibull or gamma deviate 

FPI response function 

Weibull distribution parameters 

composite thermal expansion coefficient (ppm/°F) about struc- 
tural axes 

gamma function 

gamma distribution parameters 
normal distribution parameters 
fiber and matrix Poisson's ratios 


INTRODUCTION 

The properties of the polymer matrix composites display considerable scat- 
ter because of the variation inherent in the properties of constituent materi- 
als. Distinct distributions to describe the effects of scatter on composite 
properties facilitate the composite mechanics calculations. For example com- 
posite strength is often examined probabilistically by assuming that the ply 
failure strength has a specific distribution (usually Weibull) which is then 
used in a laminate failure criterion (1). Analysis of this type has the short- 
coming that different failure mechanisms occurring at a lower level, that is, 
at the fiber and matrix level are not directly accounted for when the ply fail- 
ure stress is the primitive random variable. 

A better approach to quantify the uncertainties in the behavior of compos- 
ites would be to account for the variations in the properties starting from 
the constituent (fiber and matrix) level and integrating progressively to 
arrive at the global or composite level behavior. Typically, these uncertain- 
ties may occur at the constituent level (fiber and matrix properties), at the 
ply level (fiber volume ratio, void volume ratio, etc.) and the composite 
level (ply angle and lay-up). In this paper, a computational simulation tech- 
nique is described which accounts for uncertainties at various levels to pre- 
dict the behavior of a quasi-isotropic graphite/epoxy (0/45/90-45) s laminate. 
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MICROMECHANICAL AND MACROMECHANICAL UNCERTAINTIES 


Uncertainties at the Micromechanics Level 

To account for uncertainties at all levels of a composite, one has to 
start with uncertainties at the fiber and matrix level and use composite 
mechanics to obtain laminate level response. In the present effort the com- 
posite mechanics available in ICAN (2) is utilized to obtain the ply/laminate 
level response. At the micromechanics level 29 parameters (the constituent 
properties) are required by ICAN as input (2) (schematic Fig. 1). In addition 
three fabrication process variables (Table II) are needed to compute ply prop- 
erties. For the most part, these properties were considered to be normally 
distributed about some mean value. However, the fiber and matrix strengths 
were talcen to be distributed as a Weibull distribution which is widely accepted 
for strength distributions because of its dispersed left tail and sharp right 
tail which represents experimental data well. 

The distribution types and parameters for the fiber and matrix constitu- 
ent material properties are given in Table I. Using the Monte Carlo simula- 
tion, these distribution types reproduce histograms (frequency of occurrence) 
plots as shown in Fig. 2 for fiber longitudinal modulus and in Fig. 3 for 
fiber longitudinal strength. It would require testing of 1000 specimens to 
generate them experimentally (a rather expensive and time consuming task). It 
is worth noting that for the Weibull distribution the mean is not parameter 1 
nor is the variance parameter 2. In this case the probability density func- 
tion, mean and variances are given by (3). 


f(y) 


-Hr 1 * 



where 


X* 


1 

a 


Mean = a ^r^l + ^ , Variance = a ® 


( 1 ) 

( 2 ) 


Uncertainties at the Ply Level 

The next level of uncertainties enters at the ply level. A typical graph- 
ite fiber has a nominal diameter of 0.0003 in. which means that a single ply 
contains many fibers through the thickness. If an eight-ply graphite/epoxy 
composite is considered with a nominal thickness of 0.04 in. each ply will be 
approximately 0.005 in. thick. Taking into account an interfiber spacing of 
0.00005 (for a ply with fiber volume ratio 0.6) in. there are about 15 fibers 
through the thickness of each ply. All of these fibers will have a certain 
amount of misalignment (random orientation). To account for this randomness 
in probabilistic micromechanics, linear laminate theory is used where each ply 
is broken down (substructured) into 15 subplies (2). Each of these subplies 
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was assumed to be normally distributed about the fiber direction with fiber 
orientations lying within ±5° of the 0°-ply direction. The properties of the 
constituents are assumed to be the same in the subplies within each ply. The 
fiber volume and ply thickness were represented as normally distributed while 
the void volume was represented as a gamma distribution. A gamma distribution 
was the proper choice for the void volume ratio because there is no probability 
for zero void volume and a bias towards higher void volumes. 

As was the case for the Weibull distribution, the parameters given for 
the gamma distribution do not directly represent the mean and variance of the 
distribution. The probability density function, mean and variance for the 
gamma function are given by (3). 


ct \ X K -Xy k-1 

= flk) e y 

Mean = Variance = 

a x z 


The ply level distribution parameters are given in Table II. 


( 3 ) 

( 4 ) 


Uncertainties at the Laminate Level 

The uncertainty considered at the laminate level was that of ply orien- 
tation and thickness. Each ply in the (0/45/90/-45) s was given a normal 
distribution with a 3.33° standard deviation about the deterministic angle 
(Table III). 


MONTE CARLO SIMULATION 

Given the distributions of Tables I, II and III for fiber and matrix 
properties, ply and laminate inputs, the uncertainties in the composite 
properties need to be quantified with appropriate cumulative distribution 
functions (CDF's). One approach to achieve this is to use a Monte Carlo Simu- 
lation technique. The first step in this process involves running ICAN with 
randomly selected input variables from the predetermined probability distribu- 
tion functions many times. The output comprising of the composite properties 
is saved. The second step consists of processing the various property output 
data to compute the desired CDF. An obvious disadvantage of such an approach 
is the enormous number of output sets that must be obtained to get reasonable 
accuracy in the output CDF's. 


Generation of Normal Distributions 

In order to generate the input distributions, a uniform deviate (random 
number between 0 and 1) must first be generated. Rather than use a machine 
routine to generate the random number, a portable (machine independent) uni- 
form deviate routine from Press et al. (4) was used which was based on a three 
linear congruential generator method. This routine also had the advantage of 
being able to reinitialize the random sequence. 


4 


The uniform deviate was used to generate a normal deviate by the Box- 
Muller method (4). With the aid of the normal distribution f(y) given by 


2 

JL. 

f(y)dy - — e 2 dy 
/ Zrr 


( 5 ) 


and the transformation between uniform deviates x^,X 2 and the variables 
Y1.Y2 given by 

y x = In x 1 cos 2irx 2 , y 2 - 1 In x t sin 2 ttx 2 , (6) 


the inverse transformation can be written as 


2 2 


1 , y 2 
2ir arctan ^ 


The Jacobian of this transformation is 


3(x 1 ,x 2 ) 

3(y 1 .y 2 ) 



(7) 


(8) 


which shows that each y is distributed normally. This shows that Eq. (6) 
leads to an explicit formula for calculating a normal deviate. 


Generation of Weibull Distribution 

To generate a Weibull distribution from a uniform deviate one can inte- 
grate the probability density function and then solve for the Weibull deviate. 
This gives 


1 

y = 3 [ — In ( 1 - x)] a (9) 

as a point from the Weibull distribution where x is a uniform deviate. 


Generation of Gamma Distribution 

To generate a deviate from a gamma distribution, a uniform deviate, x, 
was taken and then the zero of the function 
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0 


x 


was found. This was numerically inefficient because it involved numerical 
integration and root finding by the bisection method, but the program ran with 
sufficient speed to overlook this fact. 

The program ICAN was modified so that the properties shown in Tables I, 

II and III were given a value from their respective distributions. Output 
for the layer and composite properties were saved for 200 samples. While 
200 samples is probably not enough to converge to the actual CDF, the results 
do show a good qualitative trend. The cumulative distribution functions were 
constructed from these samples for three typical composite properties. The 
selected composite properties are the composite longitudinal modulus E cxx , 
the composite compressive strength S cxxc , and the composite thermal expansion 
coefficient a cxx . These CDF's are shown in Figs. 4 to 6. By its symmetry, 
the CDF of E cxx appears to be normally distributed while S rYYr exhibits a 
Wei bull shape (5). 


FPI SIMULATION 

An alternative approach to obtain the required cumulative distribution 
functions is to use Fast Probability Integration (FPI) program (6). FPI helps 
generate the required CDF's quicker with reasonable accuracy and a lot less 
number of sample output data. Also, it generates more information than what 
can be expected from a Monte Carlo simulation. The additional information 
that FPI offers is the output variable sensitivity information based on the 
probabilistic inputs. 

A brief overview of FPI is given below. The reader is advised to refer 
to (6) for a detailed discussion. 

Consider a response function 

Z(X) = Z(X 1 ,X 2 ,...,X n ) (11) 

where Xi X n are random variables. Also, define the function 

g = Z(X) - Z 0 = 0 (12) 

as the limit state with Z 0 a real value of Z(X). The CDF of Z at Z 0 is 
equal to the probability that [g < 0], If the probability of a desired out- 
put, pf, is defined by 

p f = P[g<0] (13) 

an exact solution of pf can be obtained from 

Pf - ...Jf 5 (X)dX (14) 

Q “ 
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where f x (X) is the joint probability density function and Q is the region 
defined by [g < 0] . 

The evaluation of the preceding integral is often intractable and this 
leads to the need for an approximate method of evaluation p^.. In doing this, 

FPI approximates the function g using a Taylor's series expansion as a linear 

n 

*!<!!> - a o + Z) a i( U i 
i=l 

or incomplete quadratic 

n 

*2<S) ’ a o *2j a i( u i ' “i) \ 
i=l 

* 

function where is the most probable point (6) of the random variable uj . 

Note that the random variables X have been replaced by standardized normal 
variables u. The coefficients of these expansions are obtained numerically 
and then the probability P [ g< 0 ] is computed. 

Because of the approximate form of the g-function, FPI requires at least 
n+1 or 2n+l data sets to evaluate the linear or quadratic g-function coeffi- 
cients ap, ai , and bj from which the probability is found. In the present 
effort only the ply level variations in the properties (29), fiber volume ratio 
and void volume ratio are considered as random variables. This means that at 
least 32 (29 constituent properties, fiber volume ratio, void volume ratio +1) 
ICAN runs are needed for the linear approximation and 63 for the quadratic 
approximation. A typical data set to FPI consisted of ICAN run with one per- 
turbed independent variable while all others remaining at mean value. For the 
linear case, the variable was perturbed one standard deviation from its mean 
value. In the quadratic case, the independent variables were perturbed twice, 
one standard deviation each, on both sides of the mean value. 

Three typical composite properties E cxx , G CX y, and a cxx were chosen as 
the output variables for the study. Since the goal was to have a minimum 
number of ICAN runs, the CDF of E cxx was calculated using 32 data sets with 
linear FPI analysis and 125 data sets with quadratic analysis. With these two 
cases, the CDF's computed by FPI lay on top of each other indicating that 32 
data sets will give a good approximation for the CDF of E cxx . 

To identify the computational savings that FPI has over a Monte Carlo sim- 
ulation, the CDF's for E cxx , G CX y, and otg XX were compared for a 32 sample 
FPI case with a 31 and 90 sample Monte Carlo simulation (Figs. 7 to 9). It 
appears that the Monte Carlo simulation is converging to the FPI simulation, 
but the FPI simulation only needed 32 samples. 
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FPI Sensitivity Output 


As was previously mentioned, one advantage of FPI is the sensitivity data 
that is produces. Before the actual sensitivity numbers are given for the com- 
posite modulus E cxx , it will be helpful to examine how this quantity is calcu- 
lated by ICAN. 


The modulus E cxx is the (1,1) entry of the matrix (2) 


'N 


1 



S (Z li+l 


- Z 1 .)([R 1 ] T [E,][R 1 ]) 

i 


(17) 


where t c is the thickness of the composite, Zjj is the distance from the 
bottom of the composite to the ply, [Rf] is rotation matrix which is a func- 
tion of the ply angle, [E j ] is the matrix of the layer elastic constants, dis- 
tortional energy coefficient. To calculate the ply elastic moduli matrix, 

[Ef], the components are calculated from primitive variables Efu, Ef22- E ra , 
Gfi 2 , Gf 23 , Cm, kf , k m , vff 2 » and v m (2). So in the case of no ply substruc- 
turing or ply angle variation, the composite modulus should be a function of 
only these 10 primitive variables. It is noted that k m is calculated by 

kjjj = 1 — kf - ky (18) 

and is not listed as a primitive variable. Thus k v will be used instead of 
k m in the sensitivity analysis. 

For the input random variables given previously, FPI calculated a mean 
Ecxx of p « 5.744 Mpsi with a standard deviation of o = 0.363 Mpsi . The 
sensitivities at ±0.3o are given in Table IV. As would be expected, the most 
sensitive primitive variables are the fiber modulus Efu and the fiber volume 
ratio kf. Primitive variable Gf 23 has a zero sensitivity which is consistent 
with the definition of E cxx given in matrix Eq. (17). 


CONCLUSIONS 

A probabilistic evaluation of an eight-ply quasi-isotropic graph i te/epoxy 
[0/45/90/-45] s laminate was completed using two approaches. The first approach 
was to use a Monte Carlo simulation technique. The second approach was to use 
fast probability integration technique (FPI). Probabilistic inputs for this 
study included constituent micromechanical properties, fiber misalignment 
within a ply, fiber volume fraction, void volume percent and ply angle mis- 
alignment for the laminate. 

It was demonstrated that the use of the FPI program can greatly reduce 
the computations needed to generate composite CDF's. FPI was demonstrated by 
generating CDF's for E cxx , G CX y and a cxx for a graph i te/epoxy 
[0/45/90/-45] s composite in the absence of fiber misalignment. 
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The results of this investigation indicate that an integrated program 
combining ICAN and FPI is feasible. Such an integrated program offers the 
potential for a computational efficient probabilistic composite mechanics 
methodology. 
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TABLE I . - CONSTITUENT INPUT DISTRIBUTION PARAMETERS FOR I CAN 




Parameter 1 

Parameter 2 

Uni ts 

Type 

E fll 

Mpsi 

Normal 

p = 31.0 

a = 1 

.5 

E f22 

Mpsi 



p = 2 .0 

a = 0 

.10 

G f 12 

Mpsi 



p = 2.0 

a = 0 

.10 

G f23 

Mpsi 



p = 1.0 

a - 0 

.05 

v f 12 

in. /in. 



p = 0.20 

a = 0 

.01 

v f 23 

in. /in. 



p = 0.25 



a fll 

ppm/°F 



p = 0.2 



a f22 

ppm/°F 



p = 0.2 

> 


Pf 

lb/in. 3 


, 

p = 0.063 

a = 0 

003 

Nf 


Fixed 

p - 10 000 

a = 0 


df 

in. 

Normal 

p = 0.003 

a = 0 

00015 

Cf 

BTU/ lb 



p = 0.20 

o = 0.01 

Kfll 

(a) 



p - 580 

a = 2 

9 

K f22 

(a) 



p = 58 

a = 2 

9 

K f 33 

(a) 


f 

p = 58 

a = 2 

9 

SfT 

ksi 

ffeibul 1 

(3 = 400 

a = 40 

Sfc 

ksi 

Weibul 1 

13 = 400 

a = 40 

E ra 

Mpsi 

Normal 

p = 0.500 

a = 0.025 

G m 

Mpsi 

(b) 


— 


v m 

in. /in. 

Normal 

p = 0.35 

o = 0 

035 

a m 

ppm/°F 



p = 36 

a = 4 


Pm 

lb/in. 3 



p = 0.0443 

o = 0. 

0022 

Cm 

BTU/ lb 



p = 0.25 

<7=0. 

0125 

Km 

(a) 


f 

p =■ 1.25 

o = 0.06 

s mT 

ksi 

Weibul 1 

p = 15 

a = 5 


s raC 

ksi 

Weibul 1 

p = 35 

a = 20 

s mS 

ksi 

Weibul 1 

p = 13 

a =7 


Pm 

in. /in. 1% 

Normal 

■mw.i.M 

o = 0.0002 


moisture 






G m 

in . 3 /sec 

Normal 

p = 0.002 

o = 0.0001 


a = BTU • in . /hr/f t 3 /°F . 

^G m is calculated using E m and v m , and isotropy. 
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TABLE II. - PLY INPUT DISTRIBUTION PARAMETERS FOR ICAN 



Uni ts 

Distribution 

type 

Parameter 1 

Parameter 2 

kf 

Percent 

Normal 


a = 3 

k v 

Percent 

Gamma 


k = 6 

©f 

Degrees 

Normal 

mam 

a = 3.33 


TABLE III. - LAMINATE INPUT DISTRIBUTION PARAMETERS FOR ICAN 


■ 

Uni ts 

Distribution 

type 

Parameter 1 

Parameter 2 

©1 

Degrees 

Normal 

^1 D 

a = 3.33 

tl 

Inches 

Normal 

P = to 

a = 0.05t o 


TABLE IV. - NONZERO SENSITIVITY 
PARAMETERS FOR E cxx FROM 
FPI AT +0 .3a AWAY FROM ' 
MEAN OF p - 5.744 MPSI 


Primitive 

variable 

Sensi tivi ty 
parameter 

k f 

0.778 

Efll 

.624 

E f22 

.260 

Ef 12 

.130 

G m 

.060 

Em 

.036 

G f23 

.0 

v f 12 



v f 12 



v m 



k v 
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COMPONENT 



Figure 1 . - Integrated composite micro- and macromechanics 
analysis embedded in the computer code ICAN. 



Fiber Modulus, Mpsi 

Figure 2. - Monte Carlo simulation (1000 samples) of 
fiber longitudinal modulus from a normal distribution. 
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Cumulative distribution function Number of Samples 



Fiber Tensile strength, ksi 

Figure 3. - Monte Carlo simulation (1000 samples) of 
fiber longitudinal strength from a Weibull distribution. 



Figure 4. - Cumulative distribution function for composite 
([o/45/90/-45] s ) modulus Eq XX from Monte Carlo simu- 
lation (200 samples) which includes ply substructuring 
effects. 
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Compressive strength Sfc XXC ,ksi 


Figure 5. - Cumulative distribution function for composite 
([O/45/90/-45L) modulus S cxx from Monte Carlo simulation 
(200 samples) which includes ply substructuring effects. 




0 



Figure 6. - Cumulative distribution function for composite 
([o/45/90/-45] s ) thermal expansion coefficient a from 
Monte Carlo simulation (200 samples) which includes ply 
substructuring effects. 
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Cumulative distribution function 



Figure 7. - Cumulative distribution function for composite 
([o/45/90/-45] s ) modulus E cxx simulation with Monte 
Carlo (no ply substructuring) and Fast Probability 
Integration (FPI). 
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Cumulative distribution function 



Composite Shear modulas Gc Xy , Mpsi 

Figure 8. - Cumulative distribution function for composite 
([o/45/90/-45] s ) Shear modulus G cxy Simulation with 
Monte Carlo (no ply substructuring) and Fast Probability 
Integration (FPI). 
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Figure 9. - Cumulative distribution functions for the com- 
posite ([o/45/90/-45] s ) thermal expansion coefficienta CX x 
simulated with Monte Carlo (no ply substructuring) and 
Fast Probability Integration (FPI). 
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